//==============================================================================
// Entropy.cc
// ----------------------------------------------------------------------------
// Implementations of various statistical information theory computations. 
//------------------------------------------------------------------------------
// $Id: //CBT-depot/CBT++/src/Entropy.cc#1 $
//------------------------------------------------------------------------------
// $Log: $
//   
//==============================================================================
//.......1.........2.........3.........4.........5.........6.........7.........8
//

#include "Entropy.hh"

double
Entropy::FirstOrderEntropy(Sequence &S)
{
	long A=0,C=0,T=0,G=0,total = 0;

	for(int i = 0;i < S.size();i++){
		if (toupper(S[i]) == 'A') A++;
		if (toupper(S[i]) == 'C') C++;
		if (toupper(S[i]) == 'G') G++;
		if (toupper(S[i]) == 'T') T++;

		if ((toupper(S[i]) != 'N') && (toupper(S[i]) != 'X')){
			total++;
		}
	}

	double PA,PC,PT,PG;

	if (total >0){
		PA = (double)A/(double)total;
		PC = (double)C/(double)total;
		PT = (double)T/(double)total;
		PG = (double)G/(double)total;	
	}

	double h = 0;
	
	if (PA != 0) h += (PA*log(PA)/log((double)2.0));  
	if (PC != 0) h += (PC*log(PC)/log((double)2.0)); 
	if (PT != 0) h += (PT*log(PT)/log((double)2.0));  
	if (PG != 0) h += (PG*log(PG)/log((double)2.0));  
	h = -h;

	if (h == 0) h = 2;  // So, how do we want blocks on N to show up?

	return(h);
}


double
Entropy::FirstOrderEntropy(Sequence &S,int start, int end)
{
	long A=0,C=0,T=0,G=0,total = 0;

	for(int i = start;i < end;i++){
		if (toupper(S[i]) == 'A') A++;
		if (toupper(S[i]) == 'C') C++;
		if (toupper(S[i]) == 'G') G++;
		if (toupper(S[i]) == 'T') T++;
		if ((toupper(S[i]) != 'N') && (toupper(S[i]) != 'X')){
			total++;
		}
	}

	double PA,PC,PT,PG;

	if (total > 0){
		PA = (double)A/(double)total;
		PC = (double)C/(double)total;
		PT = (double)T/(double)total;
		PG = (double)G/(double)total;
	}

	double h = 0;
	if (PA != 0) h += (PA*log(PA)/log((double)2.0));  
	if (PC != 0) h += (PC*log(PC)/log((double)2.0)); 
	if (PT != 0) h += (PT*log(PT)/log((double)2.0));  
	if (PG != 0) h += (PG*log(PG)/log((double)2.0));  
	h = -h;

	if (h == 0) h=2;

	return(h);
}
